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Using fluctuating hydrodynamics we describe the slow build-up of long range spatial correlations 
in a freely evolving fluid of inelastic hard spheres. In the incompressible limit, the behavior of spatial 
velocity correlations (including r~ d -behavior) is governed by vorticity fluctuations only and agrees 
well with two-dimensional simulations up to 50 to 100 collisions per particle. The incompressibility 
assumption breaks down beyond a distance that diverges in the elastic limit. 
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In the characterization of granular matter as un un- 
usual solid, fluid or gas by Jaeger et al. this letter 
addresses the granular gas regime, controlled by inelas- 
ticity, clustering Q and collapse || . Clustering is a long 
wavelength, low frequency (hydrodynamic) phenomenon 
and inelastic collapse a short wavelength, high frequency 
(kinetic) phenomenon. In the granular gas regime, also 
called rapid granular flows, the dynamics is dominated by 
inelastic collisions. Here the methods of noncquilibrium 
statistical mechanics, molecular dynamics, kinetic theory 
and hydrodynamics are most suitable for describing the 
observed average macroscopic behavior [|J|^^,^,^|,^| and 
the fluctuations around it. 

The lack of energy conservation makes the granular 
gas, whether driven or freely evolving, behave very dif- 
ferently from molecular fluids. The essential physical pro- 
cesses and detailed dynamics are described in [|]|| and 
references therein: the similarities and differences with 
molecular fluids; lack of separation of microscales and 
macroscales, not only because the grains themselves are 
macroscopic, but also because of the existence of inter- 
mediate intrinsic scales which are controlled by the in- 
elasticity and are only well separated when the system 
is nearly elastic. A simple model which incorporates the 
inelasticity of the granular collisions, consists of inelastic 
hard spheres (IHS), taken here of unit mass and diame- 
ter, with momentum conserving dynamics. The energy 
loss in a collision is proportional to the inelasticity pa- 
rameter e — 1 — a 2 where a is the coefficient of normal 
restitution. 

For an understanding of what follows we recall two im- 
portant properties of the undriven granular gas: (i) the 
existence of a homogeneous cooling state (HCS) and (ii) 
its instability against spatial fluctuations. The hydrody- 
namic equations for an IHS-fluid, started in a uniform 
equilibrium state with temperature To, admit an HCS- 
solution (see e.g. ^,||J7)) with a homogeneous tempera- 
ture T(t), described by d t T = -2 7 owT. Here the colli- 
sion frequency is u>(T) ~ VT/lo with a mean free path Iq, 
given by the Enskog theory M for a dense system of hard 



disks or spheres (d = 2, 3) and 70 = e/2d. Then T(t) = 
To/[l + ■y uj(T )t} 2 = r exp(-27 r), where r is the av- 
erage number of collisions suffered per particle within a 
time t. It is found by integrating dr = uj{T(t))dt. More- 
over, this HCS-solution is linearly unstable once the linear 
extent L of the system exceeds some dynamic correlation 
length, which increases with decreasing e, and is propor- 
tional to h g|,H|,§||. 

The dynamics of fluctuations, say, in density, 6n(r,t), 
and flow field, u(r,i), have hardly been studied 
in sharp contrast to the large number of publications 
about the average behavior. We note that fluctuations 
are absent in hydrodynamic as well as in Boltzmann- 
Enskog-typc kinetic equations which are based on molec- 
ular chaos (mean field assumption). The objects of in- 
terest in this letter are the spatial velocity and density 
correlations 



G af) {r,t) = — J dr'(u a (r + r',t)u p (r',t)) 
G nn (r,t) = ± J dv'{5n{T + v',t)5n{v',t)), (1) 

with V — L d , and the structure factors S a ^{k.,t) and 
S nn (k, t) , which are the corresponding Fourier trans- 
forms. Goldhirsch et al. [gj initiated molecular dynamics 
studies of S nn (k, t) and S , pp (k, t) = ^aa(k, i), and re- 
lated in a qualitative way the structure at small k to the 
most unstable vorticity modes, and presented a nonlinear 
analysis to explain the enslaving of density fluctuations 
by the vorticity field Q . A more quantitative description 
of the structure factors S nn (h,t) Q and S pp (h,t) has 
been recently proposed, based on the dynamics of macro- 
scopic unstable modes (Cahn-Hilliard theory of spinodal 
decomposition ]l0|). However, numerical evidence from 
molecular dynamics for the quantitative validity of this 
theory is still lacking. 

The main goal of the present paper is to calculate the 
velocity correlation function G a p{v, t) in unforced gran- 
ular flows and to show that fluctuating hydrodynamics 
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gives a quantitative description of the spatial cor- 
relations over a large intermediate time interval, con- 
trolled by linearized hydrodynamics. Because the flows 
in freely evolving IHS systems are approximately incom- 
pressible (V • u = 0), the vorticity field, V x u, in the 
nonlinear Navier-Stokes equations is practically decou- 
pled from the other hydrodynamic fields in the system. 
This implies that the density n and the temperature 
T(t) = To cx p(— 27 t) can be considered homogeneous, 
and an approximate theory based on vorticity fluctua- 
tions alone is justified. 

Thus, we describe the Fourier modes of the vorticity 
field or transverse flow field uj_(k, t) by the mesoscopic 
Langevin equation Jul (valid for Mq < 1) 



d t u ± (k, t) + u(T(t))k 2 u ± (k, t) = f (k, t), 



(2) 



where and F are orthogonal to k, and where v{T) ~ 
IoVT is the kinematic viscosity of the IHS-fluid in the 
HCS. The random noise F is assumed to be white and 
Gaussian with a correlation 

(F a (k, t)Ff3(-k,t'))/V = B aP {T(t))k 2 S(t - t'), (3) 

and a noise strength B afi {T) = 25aaTviT)/n [[n). 
With the additional assumption (see that the IHS- 

viscosity has the same functional form as for elastic hard 
spheres or disks, the solution of the proposed Langevin 
theory provides a detailed prediction for G a p{v, t) on hy- 
drodynamic space (r > Iq) and time scales (r > 1/uj(Tq)) 
without any adjustable parameters. 

We briefly indicate how this is done by calculating 
the structure factor Sj_(k, t) = (\u± a (k, t)\ 2 )/V. Here 
the subscript a is one of the (<f — 1) equivalent trans- 
verse components of . We transform Eq. (j|) into the 
standard Langevin equation with time independent noise 
strength and coefficients. This is done by the change of 
variables u{T{t))dt = dr, u ± (k,i) = ^T(t)w(k, t) and 
F(k,f) = Lu(T{t))^/T(t){(k,T) and yields 9 T w(k,r) - 
zj_(fc)w(k, t) = f(k, t), with a growth rate z±(k) — 
7o(f — fc 2 £ 2 ) and a noise strength b a p — 2<5 Q , ( 37 £ 2 /n. 
The dynamic correlation length £ = v/lujq is time in- 
dependent and of order Zo/^/To. With the help of the 
relation (\w a (k, 0)| 2 ) = V/n, the structure factor is then 
found as 



S±(k,t) 



T(t) 



i + ex P[27o; a^)]-i ^ (4) 



which is valid for klo < 1. In the elastic limit (70 — > 0) 
the standard form of fluctuating hydrodynamics and the 
fluctuation dissipation theorem are recovered. For fc£ < 1 
this equation describes new structure of the velocity cor- 
relations G a p(Y,t) on length scales of order 27r£. At the 
end of the paper we return to the predicted structure on 
the largest scales. 



On the shortest scales (r — > 0), G a p(r,t) — > 
lT(t)/n]S a p5(r), caused by self-correlations of particles. 
As our theory describes only structure on the scale r > l Q , 
we consider the equivalent function G^~g(r, t) with the 
self-correlations substracted, which is regular at the ori- 
gin. For the same reason the structure factor S a /3(k,t) 
has a plateau value 5 a pT(t)/n for k — > 00, whereas 
S^p(k, t) — * in the same limit. The function S^ /3 (k,t) 
is an isotropic tensor field, which can be decomposed 
into two independent scalar functions of k = |k|, i.e. 
5+g(k,t) = k a kpS^(k,t) + (5 a p - k a kp)S~[(k,t), where 
a hat denotes a unit vector. The incompressibility as- 
sumption implies then Wii(k, t) = 0, and consequently 
S+(k,t) = 0. Similarly, G+j(r,t) = r a f p G+{r,t) + 

(8 a — r a rp)G~\_(r, t). Consider first the longitudinal spa- 
tial correlation Gt (r, t) = [T(t)/< d ]g||(r/£; 2 7o r), which 
is given by 



9\\{x,s) 



dq 



, ; „2 fl exp[s(l - g 2 )] - 1 

an t> - . 

I — q 2 



(5) 



where cos 9 = q • x. In the incompressible limit the 
transverse correlation function is given by g±(x,s) — 
gn(x,s) + [x/(d — l)]dg\\{x, s)/dx (see fll]], chapter 3). 
These functions can be expressed as integrals over sim- 
ple functions. As an example, we quote the result for 
d = 2: 



j' ds' exp( S ')[l - exp(-x 2 /4s')}. 



(6) 



The transverse function g±(x,s) has a negative mini- 
mum; moreover g»(x,s) is positive for all x, s,d; there 
are algebraic tails g» [x, s) ~ — (d — l)g±(x, s) ~ with 
a correction term of C(exp(— x 2 /As)). Similar algebraic 
tails occur in noneq uilibrium stationary states in driven 
diffusive systems |l2j| . These functions have structure on 
hydrodynamic space and time scales where both x = r/£ 
and s = 270T can be either large or small with respect 
to unity. At small inelasticity (70 — > 0) the dynamic cor- 
relation length and mean free path Iq are well separated. 
Details will be published elsewhere Jl3| . 

To verify the theory quantitatively we have performed 
event-driven molecular dynamics simulations of smooth 
inelastic hard disks, using square periodic boundary con- 
ditions and N — 5 x 10 3 , 2 x I0 4 and 5 x 10 4 particles. 
The inelasticity parameter e = 1 — a 2 was varied between 
0.02 and 0.8, and the area fraction from 4> = 0.02 to 0.4, 
far below the solid transition (<^ S olid = 0.665). Before 
considering the range of validity of our results, we show 
in Fig. la how the relative vorticity fluctuations w(k, r) 
have grown [SpU?]]. The HCS and linear hydrodynamics 
start to break down. The vorticity field becomes large 
and evolves into a 'dense fluid of closely packed vortex 
structures', which is still homogeneous on scales large 
compared to L v . This happens the more rapidly, the 
larger the inelasticity parameter e. 
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FIG. 1. (a) Snapshot of the momentum density at r = 80 
in a system with N = 5 x 10 4 , L = 313 (</> = 0.4), a = 0.9, 
coarse grained over cells of 6.3 x 6.3, representing a rather 
stable (still very similar at r = 160) configuration of vortex 
structures of typical diameter L v (t) ~ 25 (lo — 0.34) with 
nearest neighbors having opposite vorticity for reasons of sta- 
bility, (b) \n[E(t) / E(0)] versus r from the same simulation 
showing the linear (HCS) regime, where E(t) = T(t), with 
slope —270, up to a crossover time r cr ~ 65, and the non- 
linear cooling regime with a smaller slope, where the cooling 
through inelastic collisions is partially compensated by vis- 
cous heating. 

Apart from the restrictions to hydrodynamic space and 
time scales, there are two essential criteria limiting the 
validity of our theory: (i) System sizes L must be ther- 
modynamically large (L 3> 27r£), so that Fourier sums 
over k-space can be replaced by k-integrals. (ii) Times 
must be restricted to the linear hydrodynamic regime 
(r < r cr ), so that the system remains close to the HCS. 
Monitoring the energy per particle E(t) provides a sensi- 
tive criterion to distinguish the linear from the nonlinear 
cooling regime, where the appearance of gradients causes 
viscous heating and slows down the cooling Q (see Fig. 
lb). The crossover time r cr in Fig. lb, decreasing with 
increasing e and 4>, is an intrinsic time scale, that only 
depends on the existing gradients. At large (e, 0)-values 
criterion (ii) reduces r cr (e;0) to subhydrodynamic time 
scales, with r cr (0.5; 0.25) ~ 15 and r cr (0.3;0.4) ~ 25 as 
borderline cases. On the other hand, small (e, </>)-values 
combined with small L tend to violate criterion (i). The 
small systems with N = 5 x 10 3 , and to some extent 
even those with N — 2 x 10 4 , only satisfy the criteria (i) 
and (ii) in very narrow parameter ranges. The systems 
studied in Ref. §} {N = 1024), f§ (N = 1600) and @] 
(N = 5 x 10 3 ) are, in large regions of parameter space, 
so small that L is comparable to 2-7r£, and the periodic 
boundaries induce spurious transitions in the granular 
flows. 

We have measured in simulations the equal time spa- 
tial correlation functions G /J (r, t) with \i = {nn, ||,J-} 
by two methods, first by summing a At (vj)a /J (vj) over 
pairs of particles, where a^(vi) = {1, (v, • r), (vj • f 
and binning their relative position vectors into circular 
shells, and secondly by squaring the Fourier transform of 
the coarse-grained fields, followed by an inverse Fourier 
transformation. The second approach also provides the 



correlations in k-space, which allows us to separate the 
contributions of S\\(k,t) and S±(k,t) to Gn(r t t) and to 
G± (r, t) and test the validity of the incompressibility as- 
sumption. 
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FIG. 2. Correlation functions G M (/i = {||,_L}) versus r 
at various times, inelasticities and densities for N = 5 x 10 4 
particles, where simulation results of a single run are com- 
pared with theoretical predictions (solid line). Top row (a) 
shows 10 log[G||/T] versus 10 log r at a = 0.9 and (f> = 0.245 
(lo ~ 0.8), where r cr ~ 70. The algebraic l/r 2 -tail is clearly 
visible. Middle row (b) shows G±/T versus r sX a = 0.94 and 
4> = 0.05 (lo ~ 5.8) where r cr ~ 100. Bottom row (c) shows 
G±/T versus r at a = 0.9 and <j> = 0.4 (l ~ 0.34) where 
r cr ~ 65. Note regular oscillations with period Ro ~ 50, 
which is fixed in time. 



Figures 2a,b,c show the longitudinal and transverse 
correlation functions of the flow field of a single simu- 
lation run at N — 5 x 10 4 at several (e, (^)-values. The 
low noise levels, observed in these data and in Fig. la, 
are a consequence of the IHS collisions which have the 
tendency to make particles move parallel. There is rea- 
sonable agreement with the theoretical predictions from 
our Langevin theory, in which the viscosity is taken from 
Enskog's theory. No fitting parameters are involved. The 
longitudinal Gn (r, t) in Fig. 2a shows good agreement for 
a large range of (e, 0)-values, well beyond the linear time 
regime r cr . It exhibits the l/r 2 -tail. The minimum in 
G± (r, t) at L v (t) can be identified with the mean vortex 
diameter, and the low noise data in Fig. 2c at different r 
show that L v (t) ~ yfr is growing through vorticity dif- 
fusion. At small (e, 0)-values G±(r,t) agrees well with 
our theory, as illustrated in Fig. 2b. The l/r 2 -tail in G± 
cannot be observed in a single run because of statistical 
fluctuations. At larger densities (see Fig. 2c) one ob- 
serves small oscillations around the predicted curve with 
a characteristic length R ~ 50. The oscillations become 
more pronounced at later times, where Rq stays fixed in 
time, but varies over different runs. Comparison of G± 
at t = 80 with the snapshot in Fig. la at the same pa- 
rameters, suggests that Gj_ may be viewed as the pair 
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correlation function of a densely packed fluid of 'hard 
objects' (vortices) of typical diameter L v , the oscillation 
length i?o ~ 2L V being approximately equal to the size 
of a nearest neighbor (H — )-vortex pair. Similar complex 
structures, persisting in the nonlinear regime, are typi- 
cally observed at larger (e, 0)-values (N = 4 x 10 4 ,e = 
0.64,0 = 0.05 @; N = 5 x 10 4 ,e > 0.1,0 > 0.05; 
N = 2 x 10 4 , e > 0.05, 4> > 0.25). At smaller (e, ^-values 
(linear regime) the vortex diameter L v {t) ~ y/r, and a 
transition to a 'sheared state' is induced by the periodic 
boundaries when L v (t) ~ |Z M; for instance at r ~ 600 
for TV = 2 x 10 4 , e = 0.05, = 0.245. 

(a) (b) 
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FIG. 3. (a) Structure factors S M (fc,i) (fi = {_L, ||, nn}) 
versus k from a single simulation run at N = 5 x 10 4 , 
a = 0.9, <f) = 0.245 and r = 50 where r cr ~ 70, 
compared with the theoretical prediction (solid line) for 
S±(k,t). All structure in S/j,(k,t) contained in the interval 
fc mi „ = 2-k/L ~ 0.016 < fc < l/^ = 0.25 represents long 
range correlations of dynamic origin, (b) Correlation func- 
tions G M (r, t) (/j, = {||, _L}) from S± and (solid line), and 
from S± only (dashed line). G nn (r,t) (dotted line) corre- 
sponds tO Snn(k,t). 

The description of the velocity fluctuations G a p{v, 1) in 
this letter is based on fluctuating hydrodynamics for the 
vorticity fluctuations only, i.e. the absence of longitudinal 
fluctuations (incompressibility assumption). Fig. 3a con- 
firms that this assumption is very reasonable indeed, as 
S\\ (k, t) is vanishingly small down to very small fc-values 
(k > 1/£|| ~ 0.06). However for the smallest wavenum- 
bers, the incompressibility assumption breaks down. In 
that range the longitudinal velocity fluctuations couple 
to the second unstable mode with a dispersion re- 

lation (to second order in k) z\\(k) = 70(1 — fc 2 £|). The 
nonvanishing contributions of Su(k,t) to G a p{r,t) cause 
an exponential cut-off on length scales r > 27r£|| 3> 27r£, 
and the algebraic decay ~ 1 jr d from the vorticity mode 
represents intermediate behavior, which is well observ- 
able because the two length scales £|| and £ are in general 
quite different, e.g. £11 ~ 4.3£ in Figs. 2a,c and £11 ~ 4.5£ 
in Fig. 2b, and rapidly separate in the elastic limit. 

The results of Fig. 3a were obtained by fast Fourier 
transformation of the density and momentum fields 
(coarse-grained into 256 x 256 cells), and perform- 
ing an angular average in k-space. In the same fig- 



ure one observes that S\\(k,t) has the smallest width, 
while S±(k, t) and S nn (k,t) have a comparable width. 
Moreover, the growth rate of S±(k,t)/E(t) exceeds that 
of S\\(k,t)/ E(t), which is in turn more unstable than 
S nn (k,t). Finally, if one performs an inverse Fourier 
transform on the measured S±(k,t) and S\\(k,t) sepa- 
rately to obtain the contributions to Gu (r, t) and G± (r, t) 
(see Fig. 3b), it appears that the contributions from 
S\\ (k, t) are small, and our description of the fluctuations 
in terms of a Langevin equation based on incompressibil- 
ity is confirmed by the simulations in the linear regime 

T < T cr . 
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